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Abstract 

A new numerical method is developed for solution of the Gel'fand 
- Levitan ~ Marchenko inverse scattering integral equations. The 
method is based on the fast inversion procedure of a Toeplitz Hermi- 
tian matrix and special bordering technique. The method is highly 
competitive with the known discrete layer peeling method in speed 
and exceeds it noticeably in accuracy at high reflectance. 

1 Introduction 

Promising technological applications of fiber Bragg gratings (FBG) [1] stimu- 
late research and development of numerical methods of their synthesis. The 
propagation of counter-directional waves in single-mode fiber with quasi- 
sinusoidal refractive index modulation is described by coupled wave differ- 
ential equations [2]. Calculation of reflection coefficient r{u)) from given co- 
ordinate dependence of the refractive index is the direct scattering problem. 
The inverse scattering problem consists in recovery of the refractive index 
from given frequency dependence of the reflection coefficient r{uj). In math- 
ematical physics the inverse problem for coupled wave equations reduces to 



coupled Gel'fand - Levitan - Marchenko (GLM) integral equations [3]. How- 
ever, the straightforward numerical solution of the GLM equations is usually 
considered as too complicated for practical FBG synthesis. At first sight it 
requires A^^ operations, where is the number of discrete points along the 
grating. 

Since, the solution of integral equations seems to be inefficient, other nu- 
merical methods of FBG synthesis are elaborated. In particular, iterative 
methods with IN"^ operations are widespread, where / is the number of it- 
erations necessary for convergence. For instance, they are successive kernel 
approximations by Frangos and Jaggard [4] , high-order Born approximations 
by Peral et al [5] or advanced algorithm by Poladian [6] which uses infor- 
mation about the reflection characteristics from both ends of the grating. 
Sometimes additional approximations are applied. For example. Song and 
Shin [7] approximate the reflection spectrum by a rational function or Ah- 
mad and Razzagh [8] approximate the kernel function of integral equations 
by polynomials. 

The alternative approach is the layer peeling method known from quan- 
tum mechanics and geophysics and applied for FBG synthesis by Feced et 
al [9], Poladian [10] and Skaar et al [11]. The method has a clear physical 
interpretation of the reflected signal as a superposition of impulse responses 
from different uniform layers or point reflectors placed along the grating. 
Each thin layer has small reflectivity and can be taken into account within 
the first Born approximation. Because of high efficiency (of the order of A^^ 
operations) this method becomes widely used. The disadvantage of conven- 
tional layer peeling is the exponential decay of accuracy along the grating 
because of error accumulation during the reconstruction process [12]. The 
comparable efficiency A^^ was demonstrated by Xiao and Yashiro [13] who 
transformed the GLM integral equations to hyperbolic set of partial differ- 
ential equations and solved it numerically. This approach have several mod- 
ifications, in particular, Papachristos and Frangos [14] came to second-order 
partial differential equations and also solved them numerically. 

Better results at high reflectance are demonstrated by combination of 
the iterations and the layer peeling. It is the integral layer peeling method 
proposed by Rosenthal and Horowitz [15]. The grating is divided into thin 
layers, but layers are not assumed to have uniform proflle. The proflle of 
each layer is found by iterative solution of GLM equations. 

A recent attempt of straightforward numerical solution was made in [16]. 
The GLM equations was solved with the help of a bordering procedure and 
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Cholesky decomposition. This approach takes of the order of A^^ operations. 
The aim of present paper is to propose more efficient numerical algorithm 
with 0{N'^) operations. The improvement is possible due to specific symme- 
try of the matrix in the discrete GLM equations, the Toeplitz symmetry: the 
elements of any one diagonal are the same. The Toeplitz symmetry leads to 
considerable decrease in the number of operations, similar to the fast algo- 
rithms by Levinson [17], Trench [18] and Zohar [19]. The proposed method 
utilizes a modified bordering procedure and a second-order approximation of 
integrands, the Hermitian symmetry is also taken into account. 

The paper is organized as follows. In Sec. [2] the GLM equations are re- 
duced to convenient form for numerical calculation. The algorithm based on 
the specific "inner-bordering" technique and Toeplitz symmetry is described 
in Sec. [31 Testing numerical calculations and their comparison with the gen- 
eralized hyperbolic secant (GHS) exactly solvable profile and discrete layer 
peeling (DLP) results are summarized in Sec. HI 

2 GLM equations 

Let us consider the propagation of light through a grating with refractive 
index n + Sn{x) consisting of homogeneous background n = const and weak 
modulation 5n <C n. The refractive index modulation is quasi-sinusoidal 

5n{x)/n = 2a{x) cos {kx + 0{x)) , 

where k is the spatial frequency, a{x) is the apodization function [2] and 9{x) 
is the phase modulation that describes the chirp of the grating, variation of 
its spatial frequency. These functions are supposed to be slow-varying, that 
is, a' Ka, 6' k, where prime denotes the coordinate derivative. The 
detuning u = k — k,/2 of the light wave with respect to the grating resonance 
frequency /cq = /^/S is supposed to be small, u <^ k/2. The wave propagation 
can be described by the coupled wave equations: 

ip[ - iu^Ji = q*ip2, ip2 + 'i'^ip2 = g ^1, (1) 

where asterisk denotes the complex conjugation, the coupling coefficient q{x) 
is defined by q{x) = —ia{x)koe~^^^^\ 

The inverse problem for coupled wave equations was studied by Zakharov 
and Shabat [3], see also [7]. The problem was reduced to the Gelfand — 
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Levitan — Marchenko coupled integral equations 



Ai{x,t) + J R{t + y)A;{x,y)dy = 0, 

— oo 

X 

A2ix,t) + j R{t + y)Al{x,y)dy = -R{x + t), (2) 



X > t. 

Here 

oo 

1 



R{t) = — I r{uo)e~'^Uuj (3) 
2vr J 

— oo 

is the Fourier transform of the left reflection coefficient r{uj). For finite 
grating in the interval ^ x ^ L kernel functions Ai^2{x-,t) are not equal 
to zero only within triangular domain —x < t < x. Due to the causality 
principle the impulse respond function equals zero R{t) = at t < 0. Integral 
equations ([2]) are closed in triangular domain —x < t < x < L and allow 
one to find the kernel functions Ai^2{x, t) from function R{t) given in interval 
< t < 2L. The complex coupling coefficient q{x) can be found from the 
synthesis relation 

q{x) = 2 lim A2{x,t). (4) 

t^x—O 

For numerical analysis let us introduce more convenient variables u{x, s) = 
Al{x, X — s), v{x,t) = A2{x, T — x). GLM equations ([2]) take the form 

2x 

u{x,s) + J R*{t — s)v{x,t) dr = 0, 

s 

T 

v{x,t) + J R(t — s)u{x, s) ds = —R{t). (5) 



Functions u{x,t), v{x,t) are determined in domain ^ t ^ 2x ^ 2L. The 
synthesis relation (jl]) can be rewritten as 

q{x) =2v{x,2x -0). (6) 
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The integral operator in equations acting in the space of two-component 
vectors constructed from functions u, v is Hermitian. Note that function R in 
integrands of Eq. ([5]) depends on difference of variables only. This property 
resulting in Toeplitz symmetry of the matrix obtained by discretization of 
integral operator is exploited in the next section. 

3 Numerical procedure 

For numerical solution of Eq. ([5]) let us consider their discrete analogue. 
Divide interval ^ r ^ 2L, where function R{t) is known, by segments of 
length h = 2L/N. Introduce the discrete variables r„, s^, Xm in accordance 
with 

,m, 

, m, (7) 
,iV. 

Define grid functions u^T^ = u{xm, Tn), fn™''' = v{xm, Tn) and Rn = R{hn). 
The integrals in ([5]) can be approximated by the simplest rectangular quadra- 
ture scheme or more accurate trapezoidal scheme thus being transformed 
into sums. The accuracy of the algorithm for rectangular approximation is 
0{N~^), for trapezoidal one it is 0{N~'^). 

Discrete form of GLM equations for rectangular approximation is 

m 
n=k 

n 

V^^) +hJ2Rn-kut^ = -Rn, (8) 

k=l 

n, k = 1, . . . ,m, m = 1, . . . , N. 

The synthesis relation for the complex mode coupling coefficient ([6]) with 
accuracy 0{N~^) is 

g(™) = 2v^^\ (9) 



Sk = h { k - -] , k= 1,... 



1 
2 

mh 



T„, = hi n - - ] , n = 1, . . . 



, m = l,... 



5 



The set (jS]) at fixed index m can be represented as one matrix equation 



Q(m)^(m) ^ |g(m) 



(10) 



where vector w^™) of dimension 2m is arranged from the grid functions ui™"^ 
and vir'\ namely, 



w 



(m) 



Vector b'^'") is arranged from the zero vector of dimension m and the vector 
of dimension m with components —Rn- Square 2m x 2m matrix G*-'"-' is a 
block matrix 

G'-' ^ {^^ f ) . (11) 

Here E is the identity m x m matrix, R = R*^™) is the lower triangular 
Toeplitz m X m matrix of the form 
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Matrix R^ is the upper triangular Toeplitz mx m matrix, that is Hermitian 
conjugate to matrix R. Block matrix G'-™-' is also Toeplitz and Hermitian. 

The solution of the algebraic set fllUp can be found by the inversion of 
matrix G^*"^ using, for example, the Levinson bordering algorithm [17]. How- 
ever, we should fulfill much simpler task of finding complex mode coupling 
coefficient g*^™-* with the help of ([9]) which requires only the lower element 
of vector = to be known. Then the lower row of inverse matrix 

^Q(m)^ is interesting for us first of all. It is known that the inverse ma- 
trix to Toeplitz matrix is generally not Toeplitz, but it is persymmetric, i.e., 
symmetric with respect to the secondary diagonal [20]. Therefore, its lower 
row is the refiection of its left column 
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with respect to its secondary diagonal. The left column in its turn satisfies 
the relation 

/l\ 

Q(m)^(m) 



(13) 



The vector-column in the right hand side of (fT3ll is the first column of the 
identity matrix 2m x 2m. 

Let us also account for the Hermitian symmetry of matrix G^™-*. As 
known, the matrix inverse to Hermitian is also Hermitian. Owing to persym- 
metry and hermicity of inverse matrix its right column is 



f M 



Tilde denotes hereafter the inverted numeration of components along with 
the complex conjugation. The right column of the inverse matrix satisfies 
the relation 

/0\ 







(14) 



The last column of the identity matrix enters the right hand side. 

Since the unknown vector w{") is formed from two vectors of dimension 
m, it is convenient for us to present the left column of the inverse matrix 
f^™) as a merging of two vectors of dimension m: 



f(»») 



2;{m) 



The same relations, fll3p and f[T^ . are valid for left column f(™+i) and 
right column f of the inverse matrix (G^™"''^-*) at the next (m + l)-th 
step. 

Similar to Levinson's algorithm [17], vectors y("^+^) and z^"^^^^ at the next 
(m + l)-th step can be found by means of a bordering procedure from the 
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vectors known at the previous m-th step 



Note that the compound structure of the vectors is just what makes the bor- 
dering procedure "inner", since extending vectors y*-™'', z*^™-' by zeros means 
inserting of two rows and two columns into matrix G*-™-' with one row and 
one column placed in the middle of the matrix. At the first step we find from 
2x2 matrix G^^) that 

(1) 1 (1) — /i-Ro 
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Unknown coefficients c^, dm can be obtained from relations ([T3l) and ([1 



1 

= z — n^T^-T^' (16) 



l_|^(m)|2 ™ l_|^(m)|2 

with coefficient Z?*^™-' computed by formula 

/5M = + Rm-^yt^ + ■ • ■ + /^ii/^^) . (17) 

Then the last component v^J^^^ of vector w('"+i) is calculated as the convo- 
lution of the last row of the inverse matrix with right hand side b'^™'+^\ 

Actually, the last convolution is excessive, since relation = 2(3^^^^'^ / h 

holds. Thus, the number of arithmetic operations at each (m + l)-th step is 
of the order of m. Then the total number of required operations is A^^ which 
is approximately the same as in DLP method. 

A great advantage of the new algorithm appears when we use the trape- 
zoidal rule [21], i.e., the piecewise linear approximation of functions. The 
equations in this case remain unchanged except of the right-hand side in ([8]) 
that should be replaced by —{Rn + -Rn-i)/2 and the main diagonal of matrix 
R in f fT2|) that should be given with weight 1/2. 

Since the new procedure is based on Toeplitz symmetry of the matrix and 
the specific procedure putting a column and a row inside the matrix, we call 
it Toeplitz inner bordering (TIB) method. 
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Figure 1: The reflection spectrum of GHS grating for testing examples, koC — 
5 X 10^, JF = 3, Q = 1 (dashed hue), 2 (dots), 3 (sohd). 



4 Testing examples 

The new method is tested using a specific case of the family of exactly solvable 
chirped GHS profile of the coupling coefficient [22] 



q{x) 



(sech|)' 



It describes a FBG with apodization function 



a{x) 



sech — 

2n L 



and phase modulation 



^(a;) = 2:Fln (cosh- 



TT 



(18) 



(19) 



(20) 



where £ is the half width of grating apodization profile at level sech(l) = 
0.648, parameter Q = nCdn^a-x^/ An is the grating strength (the number of 
strokes through length L multiplied by the modulation depth of the refractive 
index). Parameter T describes the value of the chirp: the profile has a slowly 
varying spatial frequency 



de X 
Kix) = K + — = K + —T- tanh — , 
ax L L 



f211 
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Figure 2: The group delay characteristics of GHS spectrum for testing ex- 
amples at the same parameters, as in Fig. [H 



that goes smoothly from one constant spatial frequency n — IT j C to another 
n + 2T/C. 

The coupled wave equations ([H) have an exact solution that can be ex- 
pressed via the Gaussian hypergeometric function. It gives the reflection 
coefficient of the form [22] 



r[LO] 



.2-2.^2 r(^) r(/-) r(/+) 



where arguments of Euler gamma-function [23] are given by relations: 



(22) 



d 



- + i[uC- T], 



f± 



1 

i 

2 



The reflection spectrum is expressed in terms of elementary functions 
. . ,.n cosh 27r V + JF2 — cosh 27rjF 



, . (23) 

cosh 27r VQ^ + + cosh 2ttujC 

For numerical calculations we choose gratings with £ = 5 x lO^/fco, T = 3 
and Q = 1,2,3, where fco = 27™/Ao and the central resonance wavelength 
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Figure 3: Root mean square error a of the first-order (triangles) and second- 
order (boxes) reconstruction as a function of in logarithmic coordinates 
for iV = 128 4096, Q = 1, £ = 5 x WjkQ, T =\. The straight hues show 
the least-square linear fitting. 

is Ao = 1.5 /im. Their maximum reflectances, |rp = 0.6393,0.9777,0.9996, 
are referred hereafter as small, medium and high respectively. The reflec- 
tion spectrum calculated by formula (123|) is shown in Fig. [H The frequency 
detuning from resonance is shown in units 10~'^ci;o, where ioq is the central 
frequency of the reflection spectrum. The spectrum is quasi-rectangular with 
fiat top inside the Bragg reflection band. The reflectance increases with op- 
tical strength parameter Q. The width of the band Acu ^ 2V + T'^ j L 
increases, too. The group delay characteristics are plotted in Fig. [2l Each 
curve is close to straight line within the band except of the band edges. 

The GLM equations for reflection coefficient fl2^ are solved using the 
method described in Sec. [31 As the first step of calculations the fast Fourier 
transform Eq. ([3]) is performed at sufficiently long frequency interval and 
small frequency step bu) = 27i/Lmax, where Lmax = 35£, in order to neglect 
the values outside both the frequency and the coordinate intervals where the 
reflection spectrum and the grating are deflned. The frequency domain for 
integration is deflned as —Q/2 < < f2/2), where Q = N^6uj and A^^ is the 
number of discrete points in frequency. While we are going to test the method 
of solving GLM equations itself, the additional errors produced by the Fourier 
transform should be minimized. For this purpose the excessively precise 
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Figure 4: The envelope a as a function of coordinate reconstructed by second- 
order TIB (solid line) and by DLP (crosses): from the top down Q = 3,2, 1. 

determination of function R{t) is made. In order to provide the sufficient 
accuracy for the second-order method we choose iV^ 3> N, in particular, 
= 2^° at = 2^^. It does not significantly increase the total number of 
operations, since the Fourier transform requires N,^ log2 operations and 
done only once. 

The inaccuracy of rectangular and that of trapezoidal quadrature formu- 
las are compared. Root mean square error a of the grating reconstruction is 
shown in Fig. [3] as a function of 1/N. As evident from the figure the first 
and second-order algorithms result in different errors. For the first-order 
method the dependence is linear, whereas for the second order it becomes 
nearly quadratic. The slopes of fitted straight lines are 1.05 and 2.00, respec- 
tively. Moreover, the error of the second order method is significantly less at 
N ^ 2^. Then the second-order method is applied in all calculations below. 

The comparison of the second-order TIB with DLP reconstruction at 
fixed N reveals that TIB method occurs 2 4-3 times faster. The apodization 
function a{x) at the same parameters, as in Fig. [H and = 8192 is shown in 
Fig. m For relatively weak grating Q = 1,2 both methods are appropriate, 
as bottom curves demonstrate, and the resultant curves are in agreement 
with formula ([T9|) . However, for strong grating the DLP calculation gives 
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Figure 5: Comparison of the second-order TIB method with CHS profile 
(fT9|) : the deviation of numerical calculations from the analytical formula as 
a function of coordinate. The number near each curve denotes the value of 
grating strength Q. 
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Figure 6: The deviation of the spatial frequency of the grating fl?I]) from k 
calculated by TIB method (solid line) and DLP (crosses) at Q = 3. 
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significant error. The reason is probably the error amphfication in DLP [12]. 

The deviation of TIB solution from GHS profile fll9p is shown in Fig. O 
The deviation is maximal near the center of the profile and negligible at 
the ends. The curves are regular at small and medium strength and acquire 
irregular behavior for strong grating. The maximum relative error of recon- 
structed apodization function is less than 2.5-10"^ for all studied parameters. 

The phase characteristics of complex coupling coefficient demonstrate the 
similar features. At Q = 1, 2 the phase characteristics calculated by TIB and 
DLP methods are close. At high optical strength Q = 3 the error of DLP 
grows up towards the right end. The spatial frequency 9' of reconstructed 
profile is shown in Fig. O The smooth transition between two horizontal 
asymptotes of analytical expression (1^ is reproduced by TIB calculation 
for Q = 3, whereas the DLP gives the deviation at the right side of the 
curve. 

5 Discussion 

The discrete layer peeling [11] calculates q at the input end of the grating 
and then truncates the grating dealing every next step with shorter grating 
residue. This is the reason of error accumulation throughout the calculation 
from the input layer to the output one. The TIB method of matrix inversion 
recovers the complex coupling coefficient q{x) along the whole length at one 
step. Then it has higher accuracy at comparable efficiency. 

It is possible to make TIB even more efficient dividing the length L by 
segments. After reconstruction of the coupling coefficient in the current 
segment one could find the amplitudes of opposite waves at the input end 
of the next segment and repeat the procedure with the next segment. The 
efficiency could be improved if we choose the optimal number of segments. 

The similar combined procedure with indirect iterative solution of GLM 
equations, known as integral layer peeling (ILP), leads to fast reconstruction 
of a grating [15]. In that approach the grating is divided by M layers with m 
intermediate points in each. The total number of points along the grating is 

= mM. The reconstruction problem in each layer is solved by an iterative 
procedure applied to GLM integral equations. The reflection coefficient of 
truncated grating after a peeling step is found with high accuracy. The 
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computational complexity of ILP is of the order of 

ritotai ~ (^INm + log2 N 

required operations [15], where / is the number of iterations during the recon- 
struction of a layer. At I — and large m the complexity becomes less than 
N^. However with increasing m and decreasing I the accuracy goes down 
fast. 

If we were change the ILP iterations by the proposed TIB technique the 
complexity of the reconstruction within a layer would be iVln + m^. We 
obtain the total number of required operations multiplying it by M = N/m: 

"-total ~ — log2 N + mN. 

m 

This number has a minimum value minntotai ~ A^'^/^(ln A^)^/^ <^ A^^, N ^ oo 
at m ~ [NlnNY^'^ . As long as each layer is sufficiently thin and its optical 
strength is not large (Q < 1), the matrix inversion method shall give the 
superior result compared to iterations. 

For very strong gratings at 1 — |r| -^0 all the methods lose their accuracy, 
since an eigenvalue of GLM equations tends to zero and the problem becomes 
ill-conditioned. If the grating is strong, then incident light is reflected in the 
domain close to the input end. Only exponentially small part penetrates far 
from the input end, then it is almost impossible to reconstruct the profile 
of the deeper region. Fortunately, it is a formal mathematical problem. 
For more or less reasonable optical density, for instance, with maximum 
reflectance up to 99.9%, the proposed TIB method is adequately accurate. 

6 Conclusions 

Thus the new method of the FBG synthesis is proposed. The method is 
based on direct numerical solution of the coupled GLM equations. The 
Toeplitz symmetry of the matrix and the inner-bordering procedure provide 
fast computation, similar to known fast Levinson's algorithm. The second- 
order quadrature formula sufficiently improves the accuracy without loss of 
efficiency. The method is tested using exactly solvable profile of chirped grat- 
ing. The method does not concede the DLP in speed and at the same time 
remains more accurate for strong gratings. 
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